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Abstract 

The internal thermal and magnetic evolution of rocky exoplanets is critical to their 
habitability. We focus on the thermal-orbital evolution of Earth-mass planets around 
low mass M stars whose radiative habitable zone overlaps with the “tidal zone”, where 
tidal dissipation is expected to be a significant heat source in the interior. We develop a 
thermal-orbital evolution model calibrated to Earth that couples tidal dissipation, with 
a temperature-dependent Maxwell rheology, to orbital circularization and migration. 
*pedl3@uw.edu, (206) 543-0777 
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We illustrate thermal-orbital steady states where surface heat flow is balanced by tidal 
dissipation and cooling can be stalled for billions of years until circularization occurs. 
Orbital energy dissipated as tidal heat in the interior drives both inward migration 
and circularization, with a circularization time that is inversely proportional to the 
dissipation rate. We identify a peak in the internal dissipation rate as the mantle passes 
through a visco-elastic state at mantle temperatures near 1800 K. Planets orbiting a 
0.1 solar-mass star within 0.07 AU circularize before 10 Gyr, independent of initial 
eccentricity. Once circular, these planets cool monotonically and maintain dynamos 
similar to Earth. Planets forced into eccentric orbits can experience a super-cooling 
of the core and rapid core solidification, inhibiting dynamo action for planets in the 
habitable zone. We find that tidal heating is insignificant in the habitable zone around 
0.45 (or larger) solar mass stars becanse tidal dissipation is a stronger function of 
orbital distance than stellar mass, and the habitable zone is further from larger stars. 
Suppression of the planetary magnetic field exposes the atmosphere to stellar wind 
erosion and the surface to harmful radiation. In addition to weak magnetic fields, 
massive melt eruption rates and prolonged magma oceans may render eccentric planets 
in the habitable zone of low mass stars inhospitable for life. 

Keywords: tidal dissipation; thermal history; planetary interiors; magnetie field. 


1 Introduction 

Gravitational tides are common in the Solar System, from the Moon, responsible for 
driving the principle diurnal tides in Earth’s oceans and atmosphere, to lo, the most 
volcanically active body in the Solar System. Tidal dissipation as a heat source in the 
solid Earth is weak at present and often neglected from thermal history calculations 
of its interior. However, rocky exoplanets with eccentric orbits close to their star 
are expected to experience significant tides (Dole, 1964; Rasio et ah, 1996; Jackson 
et ah, 2009; Barnes et ah, 2010) that likely influence their thermal, orbital, and even 
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atmospheric evolution (Barnes et ah, 2013; Luger et ah, 2015). Recent progresses in 
modeling tidal dissipation in a visco-elastic mantle (Behounkova et ah, 2010, 2011; 
Henning et ah, 2009; Henning and Hurford, 2014) have advocated using Maxwell-type 
temperature-pressure dependent rheology and emphasized the limited applicability of 
a constant tidal quality factor “Q” model. These more complicated dissipation models 
are necessary to better characterize the tidal and orbital states of rocky exoplanets 
over a range of internal temperatures. 

The search for habitable Earth-like exoplanets commonly targets planets in orbit 
around low mass M-type stars to maximize the number of small mass planets found 
(Mayor et ah, 2014; Seager, 2013). Targeting low mass stars is beneficial for at least 
three reasons; (1) the habitable zone around M stars is much closer to the star (Koppa- 
rapu et ah, 2013), making an Earth-mass planet in the habitable zone an easier target 
for both transit and radial velocity detection, (2) low mass M stars are more abundant 
in the nearby solar neighborhood, and (3) M stars have longer main sequence times. 

On the other hand, there are several reasons why targeting M stars may be risky. 
M stars are intrinsically faint making most observations low signal to noise, and their 
flux peaks at wavelengths close to those absorbed by Earth’s atmosphere. M stars are 
more active, especially early on, which may induce massive amounts of atmospheric loss 
(Luger et ah, 2015) and biologically hazardous levels of radiation at the surface. Earth- 
mass planets in the habitable zones of M stars likely experience larger gravitational 
tides associated with star-planet and planet-planet interactions, especially considering 
that most exoplanet systems are dynamically full (Barnes and Quinn, 2004; Barnes 
and Greenberg, 2006; Van Laerhoven et ah, 2014). However, the implications of these 
tides on the thermal evolution of the interior have not yet been explored. 

The thermal history of a planet is critical to its habitability. Mantle temperature 
determines the rates of melting, degassing, and tectonics, while the thermal state of 
the core is critical to the maintenance of a planetary magnetic field that shields the 
surface from high energy radiation. The thermal evolution of the Earth and terrestrial 


3 


planets involves solving the time evolution of mantle and core temperature through a 
balance of heat sources and sinks. The thermal history of the Earth, although better 
constrained than any other planet, is still subject to significant uncertainties. However, 
avoiding both the thermal catastrophe in the mantle (Korenaga, 2006) and the new 
core paradox (Olson et ah, 2013) add significant constraints that predict a monotonic 
cooling of the mantle and an active geodynamo over the history of the planet (Driscoll 
and Bercovici, 2014). Previous models of the thermal evolution of rocky exoplanets 
(e.g. Gaidos et ah, 2010; Driscoll and Olson, 2011; Tachinami et ah, 2011; Summeren 
et ah, 2013; Zuluaga et ah, 2013) have focused on the influence of planet size on thermal 
evolution, but neglected tidal dissipation as an internal heat source, and therefore likely 
over estimate magnetic field strength and lifetime in the HZ around low mass stars. 
We improve the thermal and magnetic evolution model used in Driscoll and Bercovici 
(2014) by adding tidal heating as a internal heat source, and couple this to orbital 
evolution. 

In this paper we focus on the influence of tidal dissipation in the solid mantle of 
Earth-like exoplanets in the habitable zone around M stars. Tidal dissipation deposits 
heat in the planetary interior while simultaneously extracting energy from the orbit, 
which can lead to circularization and migration. We couple the thermal and orbital 
evolution equations into a single model to identify the conditions and timescales for 
Earth-like geophysical and magnetic evolution. Section §2 describes the thermal-orbital 
model equations. Steady state behavior is discussed in §3 to build intuition about the 
thermal-orbital coupling. Results with evolving orbits are presented in §4 and with 
fixed orbits, mimicking forcing by companion planets, in §5. The possibility of an 
internally driven runaway greenhouse is addressed in §6. The influence of tides on the 
inner edge of the habitable zone over a range of stellar masses is explored in §7. A 
summary and discussion are in §8. 
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2 Model Description 


In this section we describe the detaiis of the thermai-orbitai evoiution modei. §2.1 de¬ 
scribes the tidai dissipation modei, §2.2 describes the thermai evoiution of the coupied 
mantie-core interior, and §2.3 describes the orbitai evoiution as a function of dissipa¬ 
tion efficiency. A note on terminoiogy: here Q's refer to heat ffows (units of [W]), g’s 
refer to heat fluxes (units of [W m“^]), and the script symboi Q refers to the tidai 
quaiity factor (dimensioniess), aiso known as the specific dissipation function. 


2.1 Tidal Dissipation Model 

The gravitationai perturbation experienced by a secondary body (the pianet) in orbit 
about a primary body (the star) is approximated by the iowest order term in the 
potentiai expansion, which is the semidiurnai tide of degree 2 (Kauia, 1964). The power 
dissipated by tidal strain associated with this term in the secondary with synchronous 
rotation is (Segatz et al., 1988), 


21 , , GM^Rloje'^ 

Qtidal = —- 

2 a° 


( 1 ) 


where G is the gravitational constant, M* is stellar mass, Rp is planet radius, u is 
orbital frequency, e is orbital eccentricity, a is orbital semi-major axis, and Im(A: 2 ) is 
the imaginary part of the complex second order love number k 2 - If planetary rotation is 
synchronous then the tidal frequency is equal to the mean motion oo = n = \JGM^/a^, 
and the tidal power reduces to 

Qtidal = ~Rn{k2)G^'^Ml/^Rl^ . ( 2 ) 

This expression for tidal dissipation is the product of three physical components: (1) 
tidal efficiency (—Im(/c 2 )), (2) star-planet size (M^^'^Rp), and (3) orbit (e^/a^®/^). 

For illustrative purposes it is helpful to compare the radiative “habitable” zone 
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(HZ) to the “tidal” zone, the orbital distance at which tidal dissipation is likely to 
dominate the internal heat budget of the planet. For this comparison we compute the 
tidal heat flow using (2) for an Earth-sized planet with e = 0.1 and —Im(A: 2 ) = 3 x 10“^, 
similar to present-day Earth. Eigure 1 shows that for Mgtar < 0.3 M* the HZ overlaps 
with the tidal zone, defined as when tidal heating is the dominant term in the heat 
budget. This diagram implies that Earth-mass planets in the HZ of low mass stars 
could experience extreme tidal heating and a rapid resurfacing rate that may render 
the surface uninhabitable. The HZ and tidal zone intersect because stellar radiation 
flux is more sensitive to stellar mass than the gravitational tidal potential. We note 
that larger eccentricity or tidal dissipation efficiency (—Im(A: 2 )) would push the tidal 
zone limits out to larger orbital distances, rendering the HZ tidally dominated for larger 
mass stars. 

The one-dimensional dissipation model in (2) assumes a homogeneous body with 
uniform stiffness and viscosity. To derive the dissipation efficiency (—Im(/i; 2 )) we first 
define the love number. 


k2 


3/2 


1 + 


19 M 
2 I3st 


( 3 ) 


where /r is shear modulus and Pst is effective stiffness (Peale and Cassen, 1978). Writing 
shear modulus as a complex number and using the constitutive relation for a Maxwell 
body (Henning et ah, 2009), one can derive the dissipation efficiency in (2) as 


- Im(A:2) 
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( 4 ) 


where rj is dynamic viscosity. We note that this model does not involve a tidal Q 
factor, rather the rheological response of the mantle is described entirely by Im(/c 2 )- 
Eor comparison with other models, one can compute the standard tidal Q factor of the 
Maxwell model as 


Q = 


rjut 

Ai 


( 5 ) 
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The common approximation is then —Im(/c 2 ) ~ Re(/c 2 )/Q (e.g. Goldreich and Soter, 
1966; Jackson et ah, 2009; Barnes et ah, 2013). Although the tidal Q factor does not 
appear explicitly in the calculations below, it will be used to calibrate this model to 
Earth and is useful in comparing it to other models (also see Appendix A). 

Following previous studies (Sotin et ah, 2009; Behounkova et ah, 2010, 2011), we 
ensure that the model reproduces the observed tidal dissipation in the solid Earth by 
calibrating the effective material properties in (4) appropriately. This calibration allows 
us to approximate the total tidal dissipation over the whole mantle by a single volume 
averaged dissipation function. The Q factor of the solid Earth has been estimated 
empirically to be Qb ~ 100 (Ray and Egbert, 2012). Effective viscosity follows an 
Arrhenius Law form. 


V = Vr 


./exp 




RgTm 


/^phase 


( 6 ) 


where v = r\lpm is kinematic viscosity, pm is mantle density, v^ef is a reference viscosity, 
E,y is the viscosity activation energy, Rg is the gas constant, Tm is average mantle 
temperature, and ephase accounts for the effect of the solid to liquid phase change (see 
Table 1 for a list of constants). Shear modulus is similarly described, 





( 7 ) 


This model predicts the rapid drop in shear modulus with melt fraction demonstrated 
experimentally by Jackson et al. (2004). The reference shear modulus pref = 6.24x 10^ 
Pa and effective stiffness /3st = 1-71 x 10“^ GPa are calibrated by k 2 = 0.3 and Q = 100 
for the present-day mantle. 

We model the influence of melt fraction (j) on viscosity following the parameteriza¬ 
tion of Gosta et al. (2009), 


^phasei.^') 
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where and (jf, 'jph, and 6ph are empirical constants (Table 1). 

The functional form of tidal shear modulus at high temperature-pressure is not 
well known, so we investigate the influence of shear modulus activation energy Ap on 
the model. Contoured in Figure 2 are tidal power using (2-4), tidal power using the 
approximation —Im{k 2 ) = love number k 2 , and tidal factor Q as functions of 

shear modulus /x and viscosity v. Evolution paths of Qudai as a function of in the 
range 1500 — 2000 K are also shown in these contours for three shear modulus activation 
energies: Ap = 0, 2 x 10^, and 4 x 10® J mol~^. For the nominal activation energy 
of Ap = 2 X 10® J mol“^ the mantle cooling path passes through a maximum tidal 
dissipation around u = 10^® Pa s and fj, = 10^^ Pa, corresponding to ~ 1800 K 
where the mantle is in a visco-elastic state. In this state the short term tidal response 
of the mantle is elastic while the long-term response is viscous (see §3). We note 
that although these paths pass through a local maximum in Qudai they do not pass 
through the global maximum (dark red), which has been invoked for lo (Segatz et ah, 
1988). The main influence of increasing Ap is to shift the dissipation peak to lower 
temperatures. The nominal value of Ap = 2 x 10® J mol~^ produces a dissipation peak 
when melt fraction is about 50%. 

The tidal dissipation factor Q (Figure 2d) and tidal power using Q (Figure 2b) 
change by about one order of magnitude over the entire temperature and activation 
energy range (Figure 2d). However, tidal power using the Maxwell model in (2) fluc¬ 
tuates by 10® TW over the same range (Figure 2a). This comparison emphasizes that, 
although Q is not far from constant, dissipation using the Maxwell model is signifi¬ 
cantly different than the Q approximation (also see Appendix A). The Love number 
k 2 increases monotonically with Tm for all cases (Figure 2c) up to the limit of A ;2 = 3/2 
when ;u//3 << 2/19 in (3). 

We note that dissipation in this model is a lower bound as dissipation in the liquid 
is not included, which can occur by resonant dissipation (e.g. Tyler, 2014; Matsuyama, 
2014). Dissipation in the liquid is not likely to be a major heat source, but could drive 


mechanical flows in the core (Zimmerman et ah, 2014; Le Bars et ah, 2015) and amplify 
dynamo action there (Dwyer et ah, 2011; McWilliams, 2012). 

2.2 Thermal Evolution Model 

The thermal evolution of the interior solves the balance of heat sources and sinks in 
the mantle and core. The thermal evolution is modeled as in Driscoll and Bercovici 
(2014), with an updated mantle solidus and inclusion of latent heat release due to 
magma ocean solidification (see Appendix B.3). The conservation of energy in the 
mantle is, 


Qsurf Qconv T Qraelt Qrad “h Qcmb T Qman + Qtidal + Ql ,1 


( 10 ) 


where Qsurf is the total mantle surface heat flow (in W), Qconv is heat conducted 
through the lithospheric thermal boundary layer that is supplied by mantle convection, 
Qmeit is heat loss due to the eruption of upwelling mantle melt at the surface, Qrad is 
heat generated by radioactive decay in the mantle, Qcmb is heat lost from the core across 
the core-mantle boundary (CMB), Qman is the secular heat lost from the mantle, Qudai 
is heat generated in the mantle by tidal dissipation, and QL,man is latent heat released 
by the solidification of the mantle. Crustal heat sources have been excluded because 
they do not contribute to the mantle heat budget. Note that heat can be released from 
the mantle in two ways: via conduction through the upper mantle thermal boundary 
layer {Qconv) and by melt eruption {Qmeit)- Detailed expressions for heat flows and 
temperature profiles as functions of mantle and core properties are given in Appendix 
B. 

The thermal evolution model assumes a mobile-lid style of mantle heat loss where 
the mantle thermal boundary layers maintain Rayleigh numbers that are critical for 
convection. In contrast, a stagnant lid mantle parameterization would have a lower 
heat flow than a mobile-lid at the same temperature (e.g. Solomatov and Moresi, 2000). 
However, a stagnant lid mantle that erupts melt efficiently to the surface can lose heat 
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as efficiently as a mobile-lid mantle with no melt heat loss (Moore and Webb, 2013; 
Driscoll and Bercovici, 2014). lo is an example of this style of mantle cooling (O’Reilly 
and Davies, 1981; Moore et ah, 2007). 

Similarly, the thermal evolution of the core is governed by the conservation of energy 
in the core. 


Qcmb — Qcore Qicb “h Qrad,i 


( 11 ) 


where Qcore is core secular cooling, Qrad,core is radiogenic heat production in the core, 
and heat released by the solidification of the inner core is Qicb = Mic{Licb + Eicb), 
where Mjc is the change in inner core mass Mjc, and and Eicb are the latent and 
gravitational energy released per unit mass at the inner-core boundary (ICB). 

The internal thermal evolution equations are derived by using the secular cooling 
equation Qi = —CiMiTi, where c is specific heat and i refers to either mantle or core, 
in equations (10) and (11). Solving for Tm and Tc gives the mantle and core thermal 
evolution equations, 


Em — {Qcmb T Qrad T Qtidal T QL,man Qconv Qmelt) /-^m Cri 


( 12 ) 


Tc = -- 


{Qcmb Qrad,c) 


MqCc -AicPic^c fij' {Epe T Eq) 


(13) 


where the denominator of (13) is the sum of core specific heat and heat released by 
the inner core growth, Aic is inner core surface area, pic is inner core density, Cc is a 
constant that relates average core temperature to CMB temperature, ctRic/dTcmb is 
the rate of inner core growth as a function of CMB temperature, and Lpc and Eq are 
the latent and gravitational energy released at the ICB per unit mass (Table 1). See 
Appendix B and Driscoll and Bercovici (2014) for more details. 


2.3 Orbital Evolution Model 

The orbital evolution of the planet’s semi-major axis a and eccentricity e, assuming 
no dissipation in the primary body (the star), is (Goldreich and Soter, 1966; Jackson 
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et al., 2009; Ferraz-Mello et al., 2008) 




( 14 ) 


and 


a = 2eae. 


(15) 


Mean motion can be replaced by using = GM^^/a?, 


21 

e = —Im{k 2 ) 


M„ 


,13/2 • 


(16) 


The differential equations for thermal evolution (12, 13) and orbital evolution (16, 
15) are solved simultaneously to compute coupled thermal-orbital evolutions. 


3 Steady State Solutions 

Before exploring the full model it is useful to highlight the influence of tidal heating on 
the thermal evolution in a steady state sense by comparing heat flows as functions of 
mantle temperature. Figure 3 shows the tidal heat ffow (a) and orbital circularization 
time (b), 

^circ — ^ (^^) 

as a function of mantle temperature for a range of initial orbital distances and eccen¬ 
tricity of e = 0.1. Figure 3a shows that a peak in dissipation occurs when the mantle is 
in a partially liquid visco-elastic state (T^ ~ 1800 K), where initial tidal perturbations 
behave elastically and the long timescale relaxation is viscous. Dissipation is lower 
in a colder mantle where the response is closer to purely elastic. Dissipation is also 
lower in a hotter, mostly liquid mantle where the material behaves viscously, with little 
resistance to the external forcing. 

Also plotted in Figure 3a is the convective mantle cooling curve from (24), which 
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reflects the preferred cooling rate of the interior. Conceptually, a tidal steady state 
is achieved as the planet cools down from an initially hot state {Tm > 2000 K) until 
the convective cooling curve intersects the tidal dissipation (heat source) rate. This 
intersection implies that cooling is in balance with tidal heating so that the interior 
stops cooling. The steady state occurs around 1850 — 1950 K over the range of orbital 
distances in Figure 3a. The steady state is maintained until the orbit begins to circu¬ 
larize, which drops the dissipation curve and intersection point to lower temperatures. 
Circularization continues slowly until the dissipation rate falls below the surface cool¬ 
ing rate, at which point the planet resumes cooling normally with tidal heat playing a 
minor role in the heat budget. 

The time required to circularize, shown in Figure 3b, is inversely proportional to 
dissipation rate through (16 and 17) so that a mantle in a visco-elastic state dissipates 
orbital energy efficiently and circularizes quickly. At the inner edge of the habitable 
zone (a = 0.02 AU), circularization occurs in less than ~ 1 Gyr, while on the outer 
edge circularization requires billions of years or may not occur at all. Also shown 
in Figure 3b is the mantle cooling time It = MmCm{Tm{^) — Tm)/Qsurf, which is the 
time required for the mantle to cool from Tm{0) = 2500 K to Tm- This shows that it 
takes the mantle ~ 1 Gyr to adjust to a change in the tidal heat source. The cooling 
time is typically longer than the circularization time in the habitable zone (Figure 3b), 
implying that tidal heating can evolve faster than the thermal response of the mantle. 

4 Model Results: Evolving Orbits 

This section presents full thermal-orbital evolutions with a single Earth-mass planet in 
orbit around a 0.1 solar mass M star. We first focus on planets orbiting a 0.1 solar 
mass star, where the habitable zone is very close to the star, in order to examine an 
extreme tidal environment for a lone planet. Section 7 investigates a range of stellar 
masses. In this section the orbit of the planet is free to evolve according to (16 and 
15). Later in §5 we explore thermal evolutions with fixed orbits. The models all 
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have Afj, = 2 X 10^ J mol Tm(0) = 2400 K, and Tc(0) = 6000 K. The results are 
independent of initial mantle and core temperatures up to approximately ±500 K. 

4.1 Influence of Initial Orbital Distance 

First, we investigate the evolution of three models that start with e(0) = 0.5 at three 
orbital distances: (1) inside the inner edge of the HZ at a = 0.01 AU, (2) within the 
HZ at a = 0.02 AU, and (3) outer edge of HZ at a = 0.05 AU. 

Figure 4 compares the tidal heat flow and eccentricity of these three models as a 
function of T^. The inner case with a(0) = 0.01 AU begins with a rather high initial 
tidal heat flow {Qudai r\j 0.1 TW) considering the mantle is mostly molten at this 
time. Being so close to the star the planet has a fast circularization time (Figure 3) 
so eccentricity decreases rapidly and tidal dissipation effectively ends within 1 Myr 
(also see Figure 5c). This implies that circularization occurs during the magma ocean 
stage. 

The middle case with a(0) = 0.02 AU begins with a lower tidal heat flow because 
it orbits farther from the star. As the mantle cools and solidifies, tidal dissipation 
evolves through the visco-elastic peak at ~ 1800 K where the mantle is ~ 50% 
molten and a peak of Qudai ~ 100 TW occurs. This increase in dissipation drives a 
rapid circularization (Figure 4b), which then decreases the dissipation as the mantle 
cools further. 

The outer case at a(0) = 0.05 AU experiences the lowest initial tidal heat flow 
of Qtidai{0) ~ 10“^ TW due to it being farthest from the star. Dissipation remains 
low and the orbit remains eccentric until the mantle cools to ~ 1800 K, at which 
point dissipation increases rapidly. Tidal heat flow peaks around Qudai ~ 100 TW and 
Tm ~ 1750 K before decreasing due to decreasing eccentricity. The peak in dissipation 
occurs at a slightly lower temperature than the middle case because the eccentricity 
remains higher longer due to slower circularization. In fact, the circularization time is 
slow enough that after 10 Gyr the model still retains a finite eccentricity of e ~ 0.01 
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and a tidal heat flow of Qtidal ~ 10 TW. This shows that at the outer edge tidal 
dissipation can linger longer due to slower circularization times. 

A detailed comparison of these three models over time is shown in Figure 5. Rela¬ 
tively small differences in their temperature histories (Figure 5a) are driven by small 
differences in mantle and core cooling rates (Figure 5b). Circularization of the inner 
model occurs in the first Myr and by 100 Myr for the middle model, while the outer 
model retains a small eccentricity of e = 2 x 10“^ after 10 Gyr (Figure 5d). These 
circularization times are reflected in the tidal heat flow peaks (Figure 5b), which occur 
around 0.1 Myr for the inner case, 10 Myr for the mid case, and 1 Gyr for the outer 
case. Inward migration by 10 — 20% also accompanies this circularization (Figure 5c). 

The thermal evolutions are mainly controlled by secular cooling and radiogenic 
decay, with tidal dissipation as a temporary energy source. Mantle heating due to 
latent heat released during the solidification of the mantle is of order 10^ TW until ~ 1 
Myr, then drops below ~ 1 TW once the mantle is mostly solid around 0.1 Gyr. This 
decrease in latent heat causes the mantle heat flow to drop rapidly between 1 — 10 Myr. 
Mantle solidification, and the drop in mantle heat flow, occurs slightly later for the 
inner model for two reasons: (1) the surface heat flow is lower than the other models 
in the first Myr because the surface is hot, decreasing the upper mantle temperature 
jump; (2) tidal heating is initially moderate (~ 0.1 TW) despite the mantle being 
mostly molten due to proximity to the star. 

These mobile-lid Earth-like models have a strong temperature feedback, or ther¬ 
mostat effect, so that if mantle temperature increases (for example due to tidal dis¬ 
sipation) the viscosity decreases rapidly and the boundary layers thin out, resulting 
in an increase in the boundary heat flows. Gonsequently, increases in internal heat 
sources are accommodated by increases in heat flows so that the mantle and core cool 
monotonically. One minor exception is the brief heating of the core at 1 Gyr due to 
early radioactive decay in the core. In contrast, a stagnant lid parameterization with a 
weaker heat flow-temperature feedback would force the mantle to maintain higher tem- 
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peratures in order to accommodate the same cooling rates (e.g. Solomatov and Moresi, 
2000; Driscoll and Bercovici, 2014). Therefore, we expect mobile-lid planets to cool 
faster, dissipate tidal energy more efficiently, and circularize faster than stagnant-lid 
planets. Stagnant lid planets that are strongly tidally heated likely rely on melting 
rather than conduction to remove heat from the interior, as lo demonstrates today. 

Core cooling rates are similar at these three orbital distances, which results in 
similar magnetic moment histories and inner core nucleation times (Figure 5e). Inner 
core nucleation induces a kink in the core compositional buoyancy flux and magnetic 
moment around 4 Gyr, similar to predictions for the Earth (Driscoll and Bercovici, 
2014). Surface melt eruption rate is determined by the mantle cooling rate through 
the upper mantle geothermal gradient, so that the eruption rates at these three orbital 
distances are similar and follow the mantle heat flow history (Figure 5f). After 6 Gyr 
the middle and outer planet’s mantles are completely solid so that melt eruption ends, 
while melt eruption continues longer for the inner case due to a slightly hotter mantle. 

In summary. Earth-like planets near the inner edge of the HZ around O.IM* stars 
circularize rapidly (within a few Myr), allowing internal cooling and core dynamo ac¬ 
tion to proceed similar to Earth. On the outer edge of the HZ orbital circularization 
is slower, leading to a prolonged period of tidal dissipation that is accentuated by the 
cooling of the mantle through a visco-elastic state after ~ 1 Gyr. Despite these differ¬ 
ences in the tidal evolution, the magnetic and magmatic evolutions of these mobile-lid 
planets are similar. These three cases with high initial eccentricities of e(0) = 0.5 
demonstrate the potential for a strong coupling between orbital and thermal evolution 
by tidal dissipation. 

4.2 Summary Evolution Contours 

In this section we compare the final states (after 10 Gyr) of orbital-thermal evolution 
for 132 models that span a range of initial orbital distances of 0.01 — 0.10 AU and 
initial eccentricities of 0 — 0.5. The results are displayed as contours in initial orbital 
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a — e space (Figures 7-13). Here we consider models whose orbits evolve (left panels of 
Figures 7-13), while §5 below considers models whose orbits are fixed (right panels of 
Figures 7-13). 

Figure 6 shows the fractional change in orbital distance (a) and eccentricity (b). 
After 10 Gyr most models have circularized within the HZ due to tidal dissipation 
(Figure 6b). The iso-contour lines in Figure 6b are nearly vertical because eccentricity 
evolution is proportional to in (16). In other words, orbital circularization is 

a stronger function of orbital distance than eccentricity. Circularization also causes 
the orbit to migrate inwards (Figure 6a), although this results in a maximum inwards 
migration of only 22% of the initial distance. The evolution of orbital distance, which 
produces mainly horizontal iso-contours (Figure 6a), is controlled by initial eccentricity 
because migration is proportional to by (15); hence migration (d) is a stronger 

function of eccentricity than circularization (e). 

Figure 7a contours tidal heat flow for these models. We identify the tidal heat 
flow boundaries defined by Barnes and Heller (2013) as an Earth Twin for Qudai < 20 
TW, a Tidal Earth for 20 < Qudai < 1020 TW, and a Super-Io for Qudai > 1020 TW. 
Models that circularize by 10 Gyr have zero tidal heating. At the outer edge of the 
HZ circularization is still occurring at a rate that is proportional to the change in e in 
Figure 6b. Beyond a ~ 0.07 AU tidal dissipation is too weak to result in any signifiant 
circularization. Therefore, there are gradients in Qudai on both sides of this boundary 
at a ~ 0.07. There is also a decrease in Qudai with initial e because models with low 
initial e circularize earlier. The combination of these effects results in a peak in tidal 
dissipation around a ~ 0.07 AU and e ~ 0.5. 

This peak in tidal heat flow causes a slight increase in heat flow (Figure 8a) and 
mantle temperature (Figure 10a). This is an example of the steady state behavior 
of Figure 3 where higher tidal heat flows intersect the convective heat flow curve at 
higher mantle temperatures. Hotter mantle temperatures also cause the lower mantle 
to have lower viscosity, thinner boundary layers, and increased core heat flows (Figure 
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9a). A second maximum in core heat flow occurs at the innermost orbits due to the 
high surface temperature insulating the mantle, which keeps mantle temperature high 
(Figure 10a) and thins the lower mantle thermal boundary layer. Core temperature is 
low where core heat flow is high (Figure 11a) due to secular cooling of the core. 

After 10 Gyr of significant core cooling all models have a large solid inner core. The 
size of the inner core (also contoured in (Figure 11) is proportional to Qcmb (Figure 
9a) and is between 80 — 100% of the core radius (also see Figure 5e). Where the core 
is entirely solid no dynamo action is possible (upper left corner of Figure 12a), and 
where the core is mostly solid the magnetic moment is weak due to the small size of 
the dynamo region (upper right corner of Figure 12a). 

The eruption of melt to the surface (Figure 13a) is controlled by the upper mantle 
geothermal gradient and thus proportional to mantle heat flow with a peak around 
0.07 AU. A secondary peak in melt mass flux at close-in orbits (upper left corner of 
Figure 13a) is caused by a slightly higher mantle temperature associated with a hotter, 
insulating surface (Figure 10a). 

5 Model Results: Fixed Orbits 

In this section we consider planets whose orbits are fixed (e = d = 0). This includes 
eccentric orbits, which could be fixed, for example, through interactions with a plane¬ 
tary companion (Van Laerhoven et ah, 2014). Figures 7b-13b show contours in orbital 
space, similar to those discussed above except with fixed orbits. 

5.1 Example of Tidal Steady State 

Figure 14 shows the time evolution of a specific case with a fixed orbit of a = 0.02 
AU and e = 0.5 that reaches a tidal steady state, an example of the scenario described 
in §3. Tidal heating initially starts low (~ 10“^ TW) before increasing as the mantle 
cools for the first 10 Myr, until the mantle reaches a steady state temperature of 
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Tm ~ 1800 K. In this steady state heat loss is balanced by internal sources so that 
mantle cooling becomes insignificant. The steady state surface heat flow {Qsurf ~ 1000 
TW) corresponds to a surface heat flux similar to that of lo {qsurf ~ 2 W m“^), 
implying that this planet might be better characterized as a super-Io than Earth-like 
(e.g. Barnes et ah, 2010). The tidal steady state still allows the core to cool slowly 
because a significant temperature difference between the mantle and core persists. Core 
cooling drives a core dynamo for all 10 Gyr, although the magnetic moment rapidly 
declines as the core is nearly entirely solid by 10 Gyr (Figure 14c). 

5.2 Summary Contours 

These fixed orbit models, in contrast with the evolving models in §4, have tidal heat 
flows that are mainly determined by the orbital state rather than the cooling history. 
Specifically, Qtidal increases with e and decreases with a, producing a maximum in 
the upper left corner of Figure 7b. Mantle heat flow (Figure 8b) and temperature 
(Figure 10b) increase with tidal heating due to the positive feedback between mantle 
temperature and surface heat flow. 

Core heat flow peaks in models at moderate orbital distances where tidal heat 
flow is similar in magnitude to the sum of all other mantle heat sources; i.e. Qudai ~ 
Qcmb + Qman + Qrad (Figure 9b). This peak can be understood by considering how 
Qcmb behaves at the two tidal extremes: (1) where tidal heating is strong (upper left 
corner of Figure 7b) the mantle is forced into a hot steady state so that surface heat 
flow can accommodate all heat sources, which thins the lower mantle thermal boundary 
layer and allows a moderate core heat flow of Qcmb ~ 10 TW; (2) where tidal heating 
is weak (lower right corner of Figure 7b) the mantle and core are free to cool similar 
to Earth, so that Qcmb decreases monotonically over time. In between these limits 
tidal dissipation heats the mantle slightly, increasing the surface heat flow, but does 
not dominate the heat budget, which allows the interior to cool. Note that even 
when mantle temperature is high (~ 2000 K) it is still significantly colder than the 
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core (~ 3800 K) so that the mantle and core are not in thermal equilibrium and the 
core is forced to cool. The effect is to produce a region where a modest amount of tidal 
dissipation actually promotes core cooling, similar to the peak in the evolving models 
(Figure 9a). We refer to this 30% increase in Qcmb as the super-cooling of the core. 

The influence of this peak in core cooling rate on the dynamics of the core is 
dramatic. Intuitively, core temperature is lowest where core cooling is high (center of 
Figure 11b), and the coldest models with T^mb ~ 3850 K have completely solid cores 
(i.e. Ric = Rcore)- A fully solid core prevents fluid motion and therefore dynamo action. 
These models lose their dynamos at ~ 9.8 Gyr, so they have dynamos a vast majority 
of the time. Note that our magnetic scaling law likely provides an upper limit on 
the dynamo lifetime because the scaling law was derived from thick shell dynamos and 
does not account for stratified layers. 

This prediction implies that there is a dip (and possibly a gap) in magnetic field 
strength for tidally heated and orbitally fixed planets over most of the HZ (Figure 
12b). We note that our core liquidus does not include light element depression (e.g. 
Hirose et ah, 2013), which would tend to slow inner core nucleation and allow some 
of these models to maintain a liquid region slightly longer. This result emphasizes 
the difference between core cooling and dynamo action: cooling is ongoing (at least 
until thermal equilibration), whereas dynamo lifetime, which relies on convection in 
the liquid, is limited by the solidification time of the core (see also Gaidos et ah, 2010; 
Tachinami et ah, 2011). In other words, rapid core cooling is helpful for temporarily 
driving dynamo action but shortens the lifetime of the dynamo. 

The eruption of mantle melt to the surface follows surface heat flow (Figure 13b). 
The extreme mass fluxes of 10^® kg yr“^ correspond to a global basalt layer resurfacing 
rate of 7 m kyr“^. For reference, the Siberian traps, one of the largest igneous provinces 
on Earth and thought to be responsible for the Permian mass extinction event, is 
estimated to have produced a basalt layer at a rate of 1 m kyr“^ over the area of the 
traps (Reichow et ah, 2002). Therefore, continuous eruption rates of 10^® kg yr ^ 
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are likely to prevent such planets from being habitable. 

We also compute the same range of models but only fixing eccentricity, allowing a 
to evolve. This might occur if a neighboring planet forces the eccentricity but allows 
inward migration. In these cases we find that all models with initial orbits of a < 0.05 
AU (or a < 0.02 AU) and e > 0.2 (or e > 0.1) migrate into the central star within 10 
Gyr, and most by 5 Gyr. 

6 Internally Driven Runaway Greenhouse 

As described by Barnes et al. (2013), if interior heat flux exceeds the limit at which 
energy can be radiated from the top of the atmosphere then runaway heating of the 
surface occurs, evaporating the ocean, and leading to rapid water loss (Goldblatt and 
Watson, 2012). Figure 15 shows the time spent in an internally driven runaway green¬ 
house, defined as the period of time when the surface heat flux exceeds the threshold 
Qrunaway — 300 W m 

For both evolving (Figure 15a) and fixed (Figure 15a) orbital models the runaway 
greenhouse period is shorter at closer orbital distances, almost independent of eccen¬ 
tricity. This implies that tides, which depend strongly on eccentricity, play a minor 
role in the length of the runaway greenhouse state. The runaway greenhouse state 
is shorter for close-in planets because they have higher effective surface temperatures 
closer to the star, which insulates the mantle and decreases the initial surface heat 
flow. With lower initial surface heat flows, these inner planets drop below the runaway 
heat flow threshold earlier (Figure 5b). A second trend in Figure 15a towards even 
shorter times spent in a runaway greenhouse is found for the inner-most, high eccen¬ 
tricity planets. This drop in surface heat flow at around 50 —100 kyr occurs during the 
circularization of the inner planets’ orbits, when the tidal heat flow rapidly declines 
(Figure 5b). Gircularization causes a small dip in the surface heat flow as the interior 
temperatures and heat flows adjust to the smaller internal (tidal) heat source. This 
adjustment to lower heat flows, although seemingly minor, actually shortens the time 
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spent above the runaway threshold (Figure 5b). 

Interestingly, when the heat flow is high enough to drive a runaway greenhouse 
during the first few hundred Myr the mantle is so hot that tidal dissipation is inefficient. 
Typically tidal dissipation is not a major heat source until the mantle solidifies and 
cools down to ~1800 K, which occurs after the runaway greenhouse and magma ocean 
phases. 

In summary, we find that mobile lid Earth-like planets typically spend several 
hundred thousand years in an internally driven atmospheric runaway greenhouse state 
and that tidal dissipation in the mantle at this time plays a minor role. The runaway 
greenhouse timescale (~ 100 kyr) is shorter than the typical magma ocean solidification 
time (~ 10 Myr), a period when the surface is likely uninhabitable anyway. These 
models assume mobile lid cooling at all times, however Foley et al. (2012) proposed 
that a runaway greenhouse could induce a transition from mobile to stagnant lid, which 
would also slow internal cooling and would be detrimental to habitability. In §7 we 
explore more generally how tides may affect habitability by computing the length of 
time spent in a tidally dominated state for a range of stellar masses. 

7 Influence of Stellar Mass 

The above calculations assumed a stellar mass of 0.1 Msun- In this section we explore 
the influence of stellar mass, in the range 0.1 — 0.6 Msun, at the inner edge of the 
habitable zone. Similar to the contours in §4.2, we compute a grid of models with a 
range of initial eccentricities of 0 — 0.5 and allow the orbit to evolve in time. The initial 
orbital distance is set just outside the inner edge of the radiative habitable zone, which 
is derived from the stellar mass by the parametric equations of Kopparapu et al. (2014), 
so that the planet remains in the habitable zone after 10 Gyr of orbital migration. 

Figure 16 summarizes the results of these models in terms of two timescales: (a) 
the time spent in a tidally dominated state, dehned as when the tidal heat flow Qudai 
is 50% or more of the total surface heat flow Qsurf', (b) the time to reach Earth’s 
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present-day surface heat flow of ~ TW. 

The island-like shapes of these time contours can be explained by a combination 
of three physical effects. First, planets with initially low eccentricity (e(0) < 0.1) ex¬ 
perience weak tides and spend little time, if any, in the tidally dominated regime. At 
higher eccentricity tides become stronger, so that more eccentric planets are tidally 
dominated longer (Figure 16a). Second, as stellar mass increases the habitable zone 
moves to larger orbital distances and the tidal dissipation decreases because tidal dis¬ 
sipation in (2) is a stronger function of orbital distance (oc than stellar mass 

(oc The net result is a decrease in tidal dissipation within the habitable 

zone for increasing stellar mass, and shorter time spent in the tidally dominated state 
(Figure 16a). This effect produces contour boundaries with positive slope in Figure 
16. Third, models with high initial eccentricity (e(0) > 0.2) and close-in initial orbits 
around low mass stars {Mgtar < 0.12) experience extreme early tides that drive rapid 
orbital circularization. This leads to short times spent in the tidally dominated state. 

Figure 16b, similar to Figure 16a, shows that eccentric planets on the inner edge 
around 0.15 — 0.4 Msun stars maintain surface heat flows in excess of for 10 

Gyr due to strong tidal dissipation. Interestingly, Figure 16b shows that planets that 
experience only a temporary period of tidal heating actually cool to an Earth-like 
heat flow before 4.5 Gyr. These planets cool faster than Earth because their thermal 
adjustment timescale is longer than their circularization (or tidal heating) timescale, 
so they are still adjusting to the new heat balance with a lower tidal heat source. In 
other words, the surface heat flow, that was increased during the tidal heating phase, 
is still slightly larger than it would have been with no tidal heating. This super-cooling 
effect was also discussed in §4.2. 

In summary, tides are more influential around low mass stars. Eor example, planets 
around 0.2 Msun stars with eccentricity of 0.4 experience a tidal runaway greenhouse 
for 1 Gyr and would be tidally dominated for 10 Gyr. These time scales would increase 
if the orbits were fixed, for example by perturbations by a secondary planetary com- 
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panion. We find a threshold at a stellar mass of above which the habitable 

zone is not tidally dominated. These stars would be favorable targets in the search for 
geologically habitable Earth-like planets as they are not overwhelmed by strong tides. 

8 Discussion 

In summary, we have investigated the influence of tidal dissipation on the thermal- 
orbital evolution of Earth-like planets around M-stars with masses 0.1 — 0.6 Mgun- A 
thermal-orbital steady state is illustrated where, under certain conditions, heat from 
tidal dissipation is balanced by surface heat flow. We find that mantle temperatures 
in this balance are hotter for planets with shorter orbital distances and larger eccen¬ 
tricities. Orbital energy dissipated as tidal heat in the interior drives both inward 
migration and circularization, with a circularization time that is inversely proportional 
to the dissipation rate. The cooling of an eccentric planet in the habitable zone leads 
to a peak in the dissipation rate as the mantle passes through a visco-elastic rheology 
state. Planets around 0.1 solar mass stars with initial orbits of a < 0.07 AU circular¬ 
ize before 10 Gyr, independent of initial eccentricity. Once circular, these planet cool 
monotonically and maintain dynamos similar to Earth. Generally, we find that tidal 
dissipation plays a minor role on the dynamo history if the orbit is free to evolve in 
time. 

When the orbit is fixed the planet cools until a tidal steady state balance between 
tidal dissipation and surface cooling is reached. In the habitable zone this steady state 
can produce a super-cooling of the core when tidal heating is strong enough to heat the 
mantle and decrease its viscosity and low enough to not dominate the surface heat flow. 
This rapid cooling leads to complete core solidification, prohibiting dynamo action for 
most models in the habitable zone with e > 0.05 by 10 Gyr. In addition to weak 
magnetic fields, massive melt eruption rates in the habitable zone may render these 
fixed orbit planets uninhabitable. 

Commonly the term “habitability” refers to the influx of radiation necessary to 
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maintain surface liquid water. However, the full habitability of a planet must involve 
the dynamics of the interior and its interaction with the surface environment. We find 
that tidal heating of a planetary mantle can influence surface habitability in several 
important ways: 

1. Prolonged magma ocean stage. Close-in planets with a high eccentricity (e > 0.1) 
will experience extreme tidal heating rates of ~ 1000 TW and tidal steady state 
mantle temperatures of ~ 2000 K, implying mostly molten mantles. These super- 
tidal planets are uninhabitable as the surface itself is likely molten or close to 
the silicate solidus. 


2. Extreme volcanic eruption rates. Tidal heating, in addition to increasing sur¬ 
face heat flow, can produce extreme surface melt production rates. Even if only 
a fraction (~ 20%) of this melt erupts to the surface it can easily produce a 100 
fold increase over the present-day mid-ocean ridge eruption rate (~ 10^^ kg yr~^). 
These extreme eruption rates can lead to rapid global resurfacing and degassing 
that render the surface environment a violent and potentially toxic place for 
life. Volcanically dominated atmospheres could be significantly different from 
the modern-day Earth’s and are potentially detectable with future space- and 
ground-based telescopes (Misra et ah, 2015). 

3. Lack of magnetic field. Planetary magnetic fields are often invoked as shields nec¬ 
essary to maintain life. Magnetic fields can protect the atmosphere from stellar 
wind erosion (Driscoll and Bercovici, 2013) and the surface from harmful radia¬ 
tion (Dartnell, 2011; Griessmeier et ah, 2005). Super-cooling of the core, which 
can solidify the entire core and kill the dynamo, occurs in the habitable zone after 
~ 9 Gyr with a fixed orbit. Alternatively, a tidally heated stagnant-lid planet 
can maintain hotter mantle temperatures and lower core cooling rates, weakening 
the core generated magnetic field. Even before losing the dynamo entirely these 
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planets may have magnetic fields that are too weak to hold the stellar wind above 
the atmosphere or surface. In either case, the lack of a strong magnetic shield 
will be detrimental to life. 


4. Tidally driven runaway greenhouse. In §6 we show that driving a runaway 
greenhouse by tidal heating in the rocky interior alone is difficult. To achieve the 
runaway threshold heat flux at the surface either the planet would have to be 
forced into a highly eccentric orbit after the mantle has cooled down to ~ 1800 
K or the dissipative material properties would have to be different. For example, 
if the mantle were composed of a lower viscosity material then the maximum 
Maxwell tidal power could increase to 10^ TW (Figure 2a). A significant amount 
of tidal energy can also be dissipated in the liquid portions of the planet (Tyler, 
2014), which is beyond the scope of this study. 

With growing interest in the habitability of Earth-like exoplanets, the development 
of geophysical evolution models will be necessary to predict whether these planets have 
all the components that are conducive for life. This paper focused on a single planet 
mass, but the mathematical equations can be developed to model the evolution of other 
rocky planet/star mass ratios, including large rocky satellites around giant planets. 
However, significant uncertainties make the application to super-Earths particularly 
challenging. The fundamental physical mechanisms underpinning plate tectonics, both 
in terms of its generation and maintenance over time, are not fully understood, mak¬ 
ing extrapolation to larger planets questionable. Perhaps most importantly, material 
properties, such as viscosity, melting point, solubility, and conductivity, are poorly con¬ 
strained at pressures and temperatures more extreme than Earth’s lower mantle and 
core. This uncertainty prevails in our own Solar System where the divergent evolution 
of Earth and Venus from similar initial conditions to dramatically different present-day 
states remains elusive. 

Euture thermal-orbital modeling improvements should include coupling the evolu- 
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tion of the interior to the surface through volatile cycling and atmosphere stability. 
Advancing the orbital model to include gravitational interactions with additional plan¬ 
etary companions would allow for tidal resonances, variable rotation rates, and other 
time-dependent orbital forcings. In addition to the eccentricity tide explored here, an 
obliquity tide could also be important. Further improvements could include dissipa¬ 
tion in oceans or internal liquid layers, variable internal compositions, structures, and 
radiogenic heating rates, core light element depression, continental crust formation, 
and eventually a direct coupling of first-principles numerical simulations. 
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Figure 1: Comparison of the radiative ’’habitable” zone to the ’’tidal” zone. The radiative 
“habitable” zone is from Kopparapu et ah (2013). Inside the “tidal” zone heat released by 
tidal dissipation is likely to dominate the internal heat budget of the planet. The tidal zone 
is delineated by distances from the star where an Earth-mass planet would receive an amount 
of heat via tidal dissipation equal to either the surface heat flow of lo {Qsurf ~ 80 TW, left 
curve) or Earth {Qsurf = 40 TW, right curve). Tidal heat flow is calculated by (2) assuming 
e = 0.1 and —Im(fc 2 ) = 3 x 10“^ (^2 = 0.3, Q = 100). The gray shaded region denotes the 
zone where the planet is predicted to be radiatively “habitable” but tidally dominated, and 
therefore possibly not habitable. 
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Figure 2: Comparison of tidal properties as a function of viscosity and shear modulus for 
three shear modulus activation energies = 0 J mol~^ (blue line), 2 x 10® J moP^ (black 
line), and 4 x 10® J mol~^ (red line). Lines show tracks of z/(Tm) and //(T^) for mantle 
temperatures in the range 1500 — 2000 K. (a) Contour of tidal heat flow Qudai- (b) Contour 
of tidal power using the approximation —Im(/c 2 ) = k 2 /Q- (c) Contour of Love number ^ 2 - (d) 
Contour of tidal dissipation factor Q. Calculations use M* = O.lMsun, = 3x 10® J mol”^, 
e = 0.1. 
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Figure 3: Tidal dissipation properties as a function of mantle temperature Tm for M* = 
and e = 0.1. (a) Comparison of tidal heat flow from the Maxwell model (curves) 
with the constant Q model (Q = 100, ^2 = 0.3) at four orbital distances (see legend). 
Also shown in (a) is the mantle surface heat flow Qsurf as a function of temperature (solid 
black) and constant runaway greenhouse threshold (dashed), (b) Timescales for orbital 
circularization using the Maxwell model (same colors as in (a)) and mantle cooling (solid 
black). 
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Figure 4: Thermal-orbital evolution for three models with initial orbits of e(0) = 0.5 and 
a(0) = 0.01 (red), 0.02 (black), and 0.05 (blue). The temperature that corresponds to 50% 
melt fraction is denoted by the dashed line. 
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Figure 5: Time evolution of models with initial orbits of e = 0.5 and a = 0.01 (red), 0.02 
(black), and 0.05 (blue), (a) Temperature in the mantle (solid) and core (dashed), (b) Heat 
flow at the surface Qsurf (solid), tidal dissipation Qudai (dash-dot), mantle radiogenic heating 
(dotted), mantle latent heat (dash-dot-dot), and core heat flow Qcmb (dashed). The runaway 
greenhouse heat flow threshold (1.53 x 10^® TW) is label as a solid grey line, (c) Orbital 
distance, (d) Eccentricity, (e) Magnetic moment of core dynamo (solid) and inner core 
radius (dashed). Inner core radius axis has been scaled by core radius so the top corresponds 
to a completely solid core. For reference. Earth’s present day magnetic moment is about 80 
ZAm^. (f) Melt mass flux to the surface. Melt eruption fluxes for present-day mid-ocean 

ridges (10^^ kg yr“^) and the Siberian traps (10^^ kg yr“^) shown for reference (grey dashed). 
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Figure 6: Contour of orbital evolution after 10 Gyr for a range of initial orbital distances 
and eccentricities, (a) Change in orbital distance: (a — ao)/ao. (b) Change in eccentricity: 
(e — eo)/eo. Orbits are free to evolve in both panels. The habitable zone in denoted by 
vertical dashed white lines. 
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Figure 7: Contour of (log) tidal heat flow after 10 Gyr for a range of initial orbital distances 
and eccentricities, (a) Orbit evolves, (b) Orbit is fixed. The tidal heat flow boundaries 
defined by Barnes and Heller (2013) are shown for Earth Twins Qudai < 20 TW, Tidal 
Earths 20 < Qudai < 1020 TW, and Super-Io’s for Qudai > 1020 TW. 
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Figure 8: Contour of surface heat flow after 10 Gyr for a range of initial orbital distances 
and eccentricities, (a) Orbit evolves, (b) Orbit is fixed. White contour line shown at Earth’s 
present-day surface heat flow {Qlurf — ^0 TW). 
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Figure 9: Contour of core heat flow after 10 Gyr for a range of initial orbital distances and 
eccentricities, (a) Orbit evolves, (b) Orbit is hxed. 
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Figure 10: Contour of mantle temperature after 10 Gyr for a range of initial orbital distances 
and eccentricities, (a) Orbit evolves, (b) Orbit is fixed. 
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Figure 11; Contour of core temperature after 10 Gyr for a range of initial orbital distances 
and eccentricities. Line contours show solid core fraction Ric/Rc as a percentage (i.e. 100% 
corresponds to a completely solid core), (a) Orbit evolves, (b) Orbit is fixed. 
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Figure 12: Contour of magnetic moment after 10 Gyr for a range of initial orbital distances 
and eccentricities, (a) Orbit evolves, (b) Orbit is fixed. For reference, Earth’s present day 
magnetic moment is about 80 ZAm^. 
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Figure 13: Contour of surface melt mass flux after 10 Gyr for a range of initial orbital 
distances and eccentricities, (a) Orbit evolves, (b) Orbit is fixed. White line contour 
denotes Earth’s approximate present-day mid-ocean ridge melt flux (10^^ kg yr“^). Note 
color scales in (a) and (b) are different. 
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Figure 14: Time evolution of a model with fixed orbit of e = 0.5 and a = 0.02 AU. (a) 
Temperature in the mantle (solid) and core (dashed), (b) Heat flow at the top of the mantle 
Qconv (solid), tidal Qudai (dash-dot), mantle radiogenic heating (dotted), and core heat flow 
Qcmb (dashed), (c) Magnetic moment of core dynamo (solid) and inner core radius (dashed). 
Inner core radius axis goes from zero to total core radius. For reference. Earth’s present day 
magnetic moment is about 80 ZAm^. (d) Melt mass flux to the surface. Melt eruption fluxes 
for present-day mid-ocean ridges (10^^ kg yr“^) and the Siberian traps (10^® kg yr“^) shown 
for reference (grey dashed). 
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Figure 15: Contour of time spent in an internally driven runaway greenhouse, defined as 
when surface heat flow exceeds the threshold for a runaway greenhouse (300 W m“^, or 
1.53 X 10^^ W), for a range of initial orbital distances and eccentricities, (a) Orbit evolves, 
(b) Orbit is hxed. 
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Figure 16: Contour of (a) time spent in a tidally dominated state (i.e. Qudai/Qtotai > 0.5) 
and (b) time to reach Earth’s present-day surface heat flow {Qsurf = 40 TW). In (b) a white 
contour line is shown at 4.5 Gyr. 
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Appendix 


A Tidal Dissipation Model 

This section demonstrates the dependence of the tidal dissipation model and material 
properties on mantle temperature. Figure B.l shows several parameters related to 
the tidal dissipation rate as a function of mantle temperature for the nominal shear 
modulus activation energy of = 2 x 10® J mol~^. The Maxwell model that uses 
the full form of —Im{k 2 ) in (4) differs from the common approximation of —Im{k 2 ) ~ 
k 2 /Q for mantles hotter than the present-day (T^ > 1630 K) (Figure B.lc). The 
difference between the Maxwell model and this approximation corresponds to about 
10 orders of magnitude larger tidal heat flow at high temperature (Figure B.ld). The 
approximation is invalid at high temperature because it does not account for the drop 
in tidal dissipation expected in a liquid, since the approximation relies on Q oc rj/fi, 
which is constant, whereas the Maxwell model predicts a sharp drop in tidal dissipation 
with viscosity when n/(5 « 2/10. 

B Thermal History Model 

B.l Geotherm 

The mantle temperature profile is assumed to be adiabatic everywhere except in the 
thermal boundary layers where it is conductive. The adiabatic temperature profile in 
the well mixed region of the mantle is approximated to be linear in radius, which is a 
good approximation considering that mantle thickness D = 2891 km is much less than 
the adiabatic scale height H = Cpjag k, 12650 km, 

Tad = TuM+lad{R- r - Sum) j (18) 


37 


where the adiabatic gradient is 7 ^^ ~ 0.5 K/km. In the thermal boundary layers the 
conductive temperature solutions, 


ATj/Merf 


R — r 

^UM 


ATiMerf 


Rc - r 
Slm 


+ Ts 
Tcmb 


Upper mantle 
Lower mantle 


(19) 

( 20 ) 


replace the adiabat. Thermal boundary layer temperature jumps are ATjjm = Tjjm — 
Tg and ATlm = T^mb — Tlm, and thermal boundary layer depth is 6. Figure B.2 
shows an example whole planet geotherm T{r) at four times in the evolution. Surface 
temperature Tg is assumed to be equal to the equilibrium temperature, 

= (le^) > Pi) 


where L* is stellar luminosity and a is the Stefan-Boltzmann constant. 

The core temperature profile is assumed to be adiabatic throughout the entire 
core, i.e. the thermal boundary layers within the core are ignored. This is a good 
approximation because the low viscosity and high thermal conductivity of liquid iron 
produce very small thermal boundary layers that are insignificant on the scale of the 
whole planet. The core adiabatic profile is approximated by 

/i^2 _ ^2\ 

Tcir) = Tcmbexp —J , (22) 

where D]\f ~ 6340 km is an adiabatic length scale (Labrosse et ah, 2001). The iron 
solidus is approximated by Lindemann’s Law, 


Tpe 


Tpefi exp 




(23) 


where Tpe,o = 5600 K, 7 c is the core Gruneisen parameter, and Dpf, = 7000 km is a 
constant length scale (Labrosse et ah, 2001). This simple treatment of the core solidus 
does not account for volatile depression of the solidus, which has been demonstrated 
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experimentally (Hirose et al., 2013), and would act to slow inner core growth. Inner 
core radius can then be solved for by finding the intersection of (22) and (23). For 
details see Driscoll and Bercovici (2014). 

B.2 Mantle and Core Heat Flows 

In this section we define the remaining heat flows that appear in the mantle (10) and 
core (11) energy balance. 

The convective cooling of the mantle Qconv is proportional to the temperature 
gradient in the upper mantle thermal boundary layer, 

Qconv = ^kuM e ) (24) 

OUM 

where A is surface area and kjjM is upper mantle thermal conductivity. Qconv is written 
in terms of and the thermal boundary layer thickness 5um by requiring that the 
Rayleigh number of the boundary layer RajjM be equal to the critical Rayleigh number 
for thermal convection Rac ~ 660 (Howard, 1966; Solomatov, 1995; Sotin and Labrosse, 
1999; Driscoll and Bercovici, 2014). This constraint gives, 

Qconv = AkuM f p ^ ^ , (25) 

\RacKj 

where the thermal boundary layer temperature jump ATum has been replaced by 
~ (-UM = exp(—(i?f/M — Rm)cig/Cp) « 0.7 is the adiabatic tempera¬ 

ture decrease from the average mantle temperature to the bottom of the upper mantle 
thermal boundary layer, ATm = Tm — Tg, and the mantle cooling exponent is /3 = 1/3. 

Radiogenic heat production in the Earth is generated primarily by the decay of 
238jy^ 235fy^ 232rpj^^ 40]^^ which is approximated in the mantle by, 

Qrad{t) — Qrad,0 t/Trad) j (26) 
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where Qrad,o is the initial radiogenic heat production rate at t = 0 and Trad is the 
radioactive decay time scale that approximates the decay of the four major isotopes. 
The precise bulk silicate Earth radiogenic heat production rate is somewhat uncertain, 
so we use a nominal value of Qrad{t = 4.5 Gyr) = 13 TW (Jaupart et ah, 2007). 
Similar to the mantle convective heat flow, the CMB heat flow is, 

Qcmh = AckhM . , (27) 

olm 

where is core surface area and kiM is lower mantle thermal conductivity. The lower 
mantle and CMB temperatures, Tlm and Tcmb, are extrapolations along the mantle 
and core adiabats: Tlm = eLMTm and Tcmb = f-cTc, where eiM = exp(-(i2LM - 
Rm)oig/Cp) ~ 1.3 and Ec ~ 0.8. The lower mantle thermal boundary layer thickness 
is also derived by assuming the boundary layer Rayleigh number is critical and that 
= 2 i/c/m, which was found by Driscoll and Bercovici (2014) to produce a nominal 
Earth model. 

Core secular cooling is 

Qcore = -McCcfc , (28) 

where Me is core mass, Cc is core specific heat, and Tc is the rate of change of the 
average core temperature Tc- 

Radiogenic heat in the core is produced primarily by the decay of (Gessmann 
and Wood, 2002; Murthy et ah, 2003; Corgne et ah, 2007). Its time dependence is 
treated the same as mantle radiogenic heat in (26), but with a radioactive decay time 
scale of Trad,c = 1-2 Gyr. We assume an abundance in the core that corresponds 

to 2 TW of heat production after 4.5 Gyr. 
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B.3 Melting 

The mantle solidus is approximated by a third-order polynomial (Elkins-Tanton, 2008) , 

Tsoi{r) = Asoir^ + Bsoir'^ + Csoir + Dgoi , (29) 

where the coefficients are constants (see Table 1). This solidus is calibrated to fit the 
following constraints: solidus temperature of 1450 K at the surface, solidus temperature 
of 4150 K at the CMB (Andrault et ah, 2011), and present-day upwelling melt fraction 
of fraeit = 8%. The liquidus is assumed to be hotter by a constant offset = 500 

K, so Tiiq{r) = Tsoi{r) ATug. 

Mantle melt heat loss (or advective heat flow) is modeled as, 

Qmelt — ^erupt^melt “1“ Cm^Tmelt) , (30) 

where eempt = 0.2 is the efficiency of magma eruption to the surface (assumed to be 
constant and equal to present-day value), Mmeit is melt mass flux (see below), Lmeit 
is latent heat of the melt, Cm is specific heat of the melt, and ATmeit is the excess 
temperature of the melt at the surface (see below). This formulation of heat loss is 
similar to the ’’heat pipe” mechanism invoked for lo (O’Reilly and Davies, 1981; Moore, 
2003), where melt is a significant source of heat loss. We note that this mechanism 
is more important for stagnant lid planets where the normal conductive heat flow is 
lower (Driscoll and Bercovici, 2014). 

The melt mass flux Mmeit is the product of the upwelling solid mass flux times the 
melt mass fraction fmeit, 


Mmeit — ^upPsolidfmelti^ZUi 


(31) 


where solid density is PsoHd, volumetric upwelling rate is Vup = 1.1 QkAp/5um-, zum = 
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R — ^UM, and melt fraction is 


Mz) = ^ (32) 

liq sol 

This model predicts a ridge melt production of Mmeit = 2.4 x 10® kg s“^ for djjM = 80 
km and fmeit = 0.1, similar to present-day global melt production estimates (Cogne 
and Humler, 2004). 

We define the magma ocean as the region of the mantle with temperature exceeding 
the liquidus. Given the geotherm in (18,20) and the liquidus T)jg(r) similar to (29), 
the mantle will mainly freeze from the bottom of the convecting mantle up because the 
liquidus gradient is steeper than the adiabat (e.g. Elkins-Tanton, 2012). However, if 
the core is hot enough a second melt region exists in the lower mantle boundary layer, 
where the temperature gradient exceeds the liquidus and the mantle freezes towards 
the CMB. As can be seen in Figure B.2, a basal magma ocean exists for about 4 Gyr 
before solidifying. 

Latent heat released from the solidification of the mantle is 

QL,man — MsolLmelt , (33) 

where LmeZt is the latent heat released per kg and Mg^i is the solid mantle growth 
rate. The growth rate is calculated assuming a uniform mantle density pm so that 
^soi = PmVsoh where Vgoi = —Vug. The rate of change of the liquid volume of the 
mantle is 

Vug = , (34) 

where Tm is the mantle secular cooling rate and dVug/dTm is linearly approximated by 
8 X 10^^ m^K~^, which is the change in liquid volume from a 90% liquid to a completely 
solid mantle. This approximation implies that the latent heat released due to mantle 
solidification is linearly proportional to the mantle secular cooling rate, and the ratio 
of the latent heat flow to the mantle secular cooling heat flow is QL,manlQsec,m ~ 0.24. 
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For example, a mantle solidification time of 100 Myr corresponds to an average latent 
heat release of QL,man ~ 400 TW over that time. 

B.4 Core Dynamo 

Given the thermal cooling rate of the core, the magnetic dipole moment A4 is estimated 
from the empirical scaling law, 

M = (35) 

where 7 ^ = 0.2 is the saturation constant for fast rotating dipolar dynamos, /io = 
dvr X 10“^H m“^ is magnetic permeability, Dc = Rc — Ric is the dynamo region shell 
thickness, Rc and Ric are outer and inner core radii, respectively, and Fc is the core 
buoyancy flux (Olson and Christensen, 2006). We assume that the field is dipolar, 
ignoring the complicating influences of shell thickness and heterogeneous boundary 
conditions (e.g. Heimpel et ah, 2005; Driscoll and Olson, 2009; Aubert et ah, 2009; 
Olson et ah, 2014). In this formulation a positive buoyancy flux implies dynamo action, 
which is a reasonable approximation when the net buoyancy flux is large, but may 
overestimate the field strength at low flux. The total core buoyancy flux Fc is the sum 
of thermal and compositional buoyancy fluxes. 


Fc — Fth + F^ (36) 

where the thermal and compositional buoyancy fluxes are 


Fth 

0/-c9c ^ 

— Qc.conv 

PcCc 

(37) 

Fx 

- - 1 : (t) ■ 

(38) 


where the subscript c refers to bulk core properties, core convective heat flux is qc,conv = 
Qcmb — Qc,ad^ gravity at the ICB is approximated by gic = gcRic/Rc, and the outer core 
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compositional density difference is = pc — Px with the light element density. 
For simplicity, the expression for light element buoyancy (38) ignores buoyancy due to 
latent heat release at the ICB because it is a factor of 3.8 less than buoyancy of the 
light elements. 

The isentropic core heat flux at the CMB, proportional to the gradient of (22), is 

QcM = kcTcmbRc/D% , (39) 

where core thermal conductivity is approximated by the Wiedemann-Franz law, 

kc — O'cLcTcmb ) (^^) 

and electrical conductivity is ac and Lc is the Lorentz number. For typical values of 
high pressure-temperature iron, cjc = lOx 10^ m“^ (Pozzo et ah, 2012; Gomi et ah, 

2013), Lc = 2.5 X 10“® WnK“^, and Tcmb = 4000 K, the core thermal conductivity is 
kc = 100 Wm-^K-^ 
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Figure B.l: Tidal dissipation parameters as a function of mantle temperature for the 
nominal shear modulus activation energy oi = 2x 10® J mol~^. Circles denote the chosen 
calibration points for the present-day mantle. 
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Figure B.2: Geotherm temperature profile T(r) for the orbit evolving model with initial 
orbit of a = 0.05 AU and e = 0.5. Temperature profiles are shown at 1 Myr, 1 Gyr, 4.5 Gyr, 
and 10 Gyr. The core liquidus and mantle liquidus and solidus are shown in red. 
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Symbol 

Value 

Units 

Reference 


3 X 10^ 

J mol” 3 

Viscosity activation energy in (6) 


2 X 10^ 

J mol” 3 

Nominal shear modulus activation energy in (7) 

^sol 

-1.160 X 

CO 

Solidus coefficient in (29) (ET08) 

a 

3 X 10“® 

K-i 

Thermal expansivity of mantle 

ttc 

1 X 10-5 

K-i 

Thermal expansivity of core 

B 

2.5 

nd 

Melt fraction coefficient in (8) 

Bgoi 

1.708 X 10"^ 

K/m^ 

Solidus coefficient in (29), calibrated 

13 

1/3 

nd 

Convective cooling exponent in (25) 

(3st 

1.71 X 10^ 

GPa 

Effective mantle stiffness, calibrated in §2.1 

Cm 

1265 

J kg-i K-3 

Specific heat of mantle 

C-c 

840 

J kg-i K-i 

Specific heat of core 

Csol 

-9.074 X 10-3 

K/m 

Solidus coefficient in (29), calibrated 

D 

2891 

km 

Mantle depth 

Dpe 

7000 

km 

Iron solidus length scale 

Dn 

6340 

km 

Core adiabatic length scale 

13 sol 

1.993 X 10^ 

K 

Solidus coefficient in (29), calibrated 

3ph 

6 

nd 

Rheology phase coefficient in (8, 9) 

Eg 

3 X 10^ 

J kg-i 

Gravitational energy density release at ICB 

euM 

0.7 

nd 

Upper mantle adiabatic temperature drop 

^LM 

1.3 

nd 

Lower mantle adiabatic temperature jump 

ec 

0.8 

nd 

Average core to CMB adiabatic temperature drop 


0.8 

nd 

Rheology phase coefficient in (8, 9) 

guM 

9.8 

_9 

ms 

Upper mantle gravity 

qlm 

10.5 

_9 

ms 

Lower mantle gravity 

9c 

10.5 

_9 

ms 

CMB gravity 

Ic 

1.3 

nd 

Core Gruneisen parameter 

'Idip 

0.2 

nd 

Magnetic dipole intensity coefficient in (35) 
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Symbol 

Value 

Units 

Reference 

Iph 

6 

nd 

Rheology phase coefficient in (8, 9) 

kuM 

4.2 

W m-i K-i 

Upper mantle thermal conductivity 

kiM 

10 

W m-3 K-i 

Lower mantle thermal conductivity 

K 

10® 

2 —1 
m s 

Mantle thermal diffusivity 

Lpe 

750 

kJ kg“^ 

Latent heat of inner core crystallization 

Lmelt 

320 

kJ kg“^ 

Latent heat of mantle melting 

Le 

2.5 X 10"® 

w n K-i 

Lorentz number 

L=^ 

3.09 X 10^3 

W 

Stellar luminosity for M* = Q.lMgun (B13) 

Mm 

4.06 X 10^4 

kg 

Mantle mass 

Me 

1.95 X 1024 

kg 

Core mass 

IJ-ref 

6.24 X 10^ 

Pa 

Reference shear modulus in (7) 

1^0 

Ap X 10“^ 

H m“^ 

Magnetic permeability 

^ref 

6 X 10^ 

2 —1 
ms 

Reference viscosity 


2 

nd 

Viscosity jump from upper to lower mantle 

Qrad,0 

60 

TW 

Initial mantle radiogenic heat flow (J07) 

R 

6371 

km 

Surface radius 

Rc 

3480 

km 

Core radius 

Rm 

4925 

km 

Radius to average mantle temperature Tm 

Rdc 

660 

nd 

Critical Rayleigh number 

pc 

11900 

kg m“3 

Core density 

pic 

13000 

kg m“3 

Inner core density 

Pm 

4800 

kg m“3 

Mantle density 

Pmelt 

2700 

kg m“3 

Mantle melt density 

Psolid 

3300 

kg m“3 

Mantle upwelling solid density 

^Px 

700 

kg m“3 

Outer core compositional density difference 

CTc 

10 X 10® 

S m“^ 

Core electrical conductivity 

Tpefl 

5600 

K 

Iron solidus coefficient in (23) 
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Symbol 

Value 

Units 

Reference 

'^rad 

2.94 

Gyr 

Mantle radioactive decay time scale 

'^rad^c 

1.2 

Gyr 

Core radioactive decay time scale 


5 X 10"^ 

nd 

Rheology phase coefficient in (8, 9) 


Table 1: Model constants. Non-dimensional nnits are denoted n.d. References 
are: B13=Barnes et al. (2013); ET08=Elkins-Tanton (2008); H00=Hirsclimann 
(2000); J07=Janpart et al. (2007); M84=McKenzie (1984); M88=McKenzie and 
Bickle (1988). 
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